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Abstract 

We study an example of exact parametric resonance in a extended system 

ruled by nonlinear partial differential equations of nonlinear Schrodinger type. 

It is also conjectured how related models not exactly solvable should behave 

in the same way. The results have applicability in recent experiments in 

Bose-Einstein condensation and to classical problems in Nonlinear Optics. 
PACS: 03.40.Kf, 03.75.-b 42.65.-k, 
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Resonances are one of the recurrent themes of Physics. After reading the current PACS 
sectioning scheme one finds that the word resonance appears in fifty different items which 
is only a naive way to measure one very important concept in Physics. 

When speaking of resonances we use to refer to an "anomalously large response to a 
(maybe small) external perturbation". The best known model is the undamped harmonic 
oscillator driven by an external force, which is textbook material. In this trivial system one 
easily proves existence of unbounded linearly increasing terms in the solution. A more com- 
plicated kind of resonant behavior are parametric resonances, where a relevant parameters of 
the system is modulated to achieve resonance. This phenomenon has been known since the 
Middle Ages and its simplest examples are the parametrically forced harmonic oscillator 
modeled by Hill's equation x + p(t)x = [0] and a particular version, p(t) = 1 + e cos cut, 
known as Mathieu's equation. 

Since those are linear problems much can be said on their solutions. For instance, 
Mathieu's equation can be studied by means of Floquet's theory for periodic coefficient 
equations and resonance existence can be proven analytically -even in the presence of linear 
dissipation- for various parameter regions on the (e, u) plane. However resonance phenomena 
are not restricted to linear problems but also appear in nonlinear finite dimensional problems, 
e.g. in many Hamiltonian chaotic systems where torus resonances are the reason for the 
origin of chaos, in impact oscillators and in nuclear magnetic and spin resonances, to cite 
only a few examples. There is finally the framework of systems modeled by partial differential 
equations. Here some results are known in simple linear cases, but in general the study of 
resonances in extended (i.e. ruled by partial differential equations) nonlinear systems poses 
many open questions. 

Our aim in this letter is to prove that unbounded resonances are possible in simple, 
experimentally realizable Hamiltonian wave problems, specifically in the framework of non- 
linear Schrodinger (NLS) equations [Q]. This is a very surprising result since one is tempted 
to think (as it usually happens) that a resonant perturbation can only excite a small num- 
ber of modes of the infinite present in an extended problem. Since the nonlinearities mix 
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different modes the amplification of one of them is stopped by energy transfer to the oth- 
ers, specially in conservative systems where the conservation laws control the number and 
amplitude of active modes. Thus, it could be thought that once the energy is transfered to 
the non-resonant modes and being their growth limited by nonlinear constraints, the action 
of the perturbation should be controlled and unbounded growth of the relevant quantities 
should be inhibited. Indeed this is the common mechanism in many systems. We will show 
that it is not the only possibility and truly resonant behavior is possible even in simple 
physically relevant models. 

The model.- Let us consider the following Nonlinear Schrodinger equation in n spatial 
dimensions. 

i^- = -^Au+\u\ 2 u + V(r,t)u, (1) 

where A = Y^=i d 2 /dx 2 . This equation is the adimensional description of many different 
phenomena |3|. One of the hottest topics where it appears is in the mean field model 
of nonuniform trapped Bose-Einstein condensates (BEC) j|] where V(r) = | J21=i ^j(t)tf- 
Time dependent coefficients reflect experiments in which the trap is perturbed to obtain the 
response spectrum of the condensate []5| (For a theoretical analysis of these experiments see 
Ref. H). A more detailed study of the condensate response to time dependent perturbations 
was performed in f?J and the existence of strong resonances was proposed on the ground of 
approximate variational methods and numerical simulations of the Eq. (|l]). Nonetheless, the 
nature of the resonance could not be completely assured on the ground of approximate or 
numerical techniques since both have limited applicability, specially for the strong response 
phenomena considered here. 

Let us first consider in this letter the two-dimensional case with time-dependent parabolic 
potential 

= ~ Au + \u\ 2 u + ±\(t) (x 2 + y 2 ) u, (2) 

where A = d 2 /dx 2 + d 2 /dy 2 . Although Bose-Einstein condensation problems are usually 
three dimensional, this model can be related to pancake type traps, and to recent experiments 
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with quasi two-dimensional condensates or coherent atomic systems The model (||) is 

also known in nonlinear Optics where it is used to study the propagation of paraxial beams 
in fibers with a (modulated) parabolic profile of the refraction index. 

From the mathematical point of view, Eq. (0) is a parametrically forced NLS equation 
that has not been studied rigorously before. Related systems are the externally driven 
damped NLS in one dimension ||10|| , the AC-driven damped sine-Gordon system and one 
type of non-self-adjoint parametrically driven NLS [pL2|,|T3l| . However those works concentrate 
on one dimensional equations and mostly on soliton stability problems. 

Moment method.- We will first study the radially symmetric version of Eq. (0), search- 
ing solutions of the form ip(r,6,t) = u(r,t)e tme , which includes both the typical radially 
symmetric problem corresponding to m = 0, and vortex line solutions, with m ^ 0. The 
simplified equation for u is 



.du 1 d / du\ + ^m 2 + ^ ^ ^ \{t) 
2r dr \ dr J 



r 2 u, 



(3) 



dt 2r dr \ dr J \ 2r ' 1 ~ 1 ' 2 

Eq. (^|) is non-integrable and has no exact solutions even in the constant A case. Its solutions 
are stationary points of the action S = J/ 1 £(t), where 



u- 



du* 
~dt 



1 

+ 2 



du* 



Or 



*9u\ , 2 
u — \ a x 
dt 



m o 
H — -\u\ + \u 



X(t) 



2 

' 1 d 2 x. 



r \u\ a x 



(4) 



Let us define the following integral quantities [14 
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where d 2 x = 2irrdr because of the symmetry. These magnitudes are related physically to 
the norm (intensity or number of particles), width, radial momentum and energy of the wave 



packet. In the Optical interpretation of Eq. (Bl) these quantities are known as moments and 



are used in (usually approximate) calculations related to beam parameters evolution flH,l5 
It is remarkable that the Ij satisfy simple, and most important, closed evolution laws 

f - <»> 

dli 

^ = -2A(i)/ 2 + 4/ 4 , (6c) 

f = ~\^- ™ 
The first equation comes from the phase invariance of Eq. (^) under global phase transforma- 
tions and corresponds to the L 2 -norm conservation (in BEC it is interpreted as the particle 
number conservation in the mean field model). The other equations can also be obtained in 
connection with the invariance of the action under symmetry transformations. 

Singular Hill's equation.- Eqs. @ form a linear non-autonomous system for the un- 
knowns, Ij,j = 1,...,4 that has several positive invariants under time evolution, of which 
the most important is 

Q = 21 J 2 - if /4 > 0. (7) 

With the help of this quantity, the system (|6|) can be reduced to a single equation for the 
most relevant parameter I^if), which is 

If we were able to solve Eq. then the use of Eqs. (^) would allow us to track the evolution 
of the other ones. We can do so by defining X(t) = V^2? whose physical meaning is the 
wavepacket width, and substituting it into (H). This procedure gives us 

X + X(t)X = ^- r (9) 

The resulting equation is a singular (nonlinear) Hill equation. In general it cannot be solved 
explicitly (see Ref. [RJ for other examples of these equations), however in this case we are 
fortunate that the general solution of Eq. (^]) can be obtained [|l^] and is given by 



m = sfm + ^(t), (io) 

where u(t) and v(t) are the two linearly independent solutions of the equation 

x + X{t)x = 0, (11) 

which satisfy u(to) = X(to),ii(t ) = X'(to),v(t ) = 0,v'(t ) ^ 0, and W is the Wronskian 
W = uv — uv = const. ^ 0. 

The conclusion is that the evolution of the width of the wave packet, which is given by 
( PD is closely related to the solutions of the Hill equation (PP. This is a very well studied 
problem || which is explicitly solvable only for particular choices of X(t), but whose solutions 
are well characterized and many of its properties are known. 

A practical application of Eq. ( |10"D is that one may design X(t) starting from the desired 
properties for the wave packet width evolution. This can be done in BEC applications by 
controlling the trapping potential and in Optics by the precise control of the z— dependent 
refraction index of the waveguide. 

If it is supposed that X(t) depends on a parameter X(t) = 1 + X(t) with X(t) a periodic 
function with zero mean value and peak value e (not necessarily small), then there is a 
complete theory describing intervals of e for which all solutions of Eq. ( pTTj ) are bounded 
(stability intervals) and intervals for which all solutions are unbounded (instability intervals). 
Both types of intervals are ordered in a natural way. Let us finally concentrate in the 
physically relevant case X(t) = 1 + e cos ujt, and where e needs not to be a perturbative 
parameter. In that case the parameter regions where exact resonances occur can be studied 
by several means. First, for any fixed e, there exist infinite ordered sequences {uj n }, {u>' n } 
tending to zero such that Eq. (|TT[) is resonant if uj belongs to (u n ,u' n ) for some n. Second, 
when uj is fixed, resonances appear for e big enough. Further, a stability diagram can be 
drawn in the e — uj plane as shown in Fig. |]. The boundaries of these regions are the so- 
called characteristic curves, which can not be explicitly derived but whose existence can be 
proven analytically: if D(e,uj) is the discriminant of the equation, characteristic curves are 
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obtained by solving the equations D(e,u) = 2 and D(e,u) = —2. In particular, instability 
regions start on the frequencies u; = 2,l,l/2,...,2/n 2 ,... Q [Fig. |]. 

The resonant behavior depends only on the mutual relation of the parameters but not on 
the initial data. Since a resonance in Eq. ( [TT]) automatically implies unbounded behavior 
in the solution [Eq. flT0|)] we obtain the result that for resonant parameters all solutions 
have divergent width provided they start from finite values of the moments (which is the 
generic case). With respect to stability, it is deduced from Massera's Theorem that if (e,u) 
belongs to a stability region of Eq. (0), there exists a periodic solution of Eq. (J|) which is 
Lyapunov stable. 

We emphasize that our analysis of the cilindrically symmetric problem performed up to 
now is completely rigorous. Now we will turn to approximations to discuss other related 
problems. 

N on- symmetric problems.- Throughout this paper we have made use of a special sym- 
metry for the confining potential and for the solution in order to make our theory exact. 
When these constrains are removed and we turn to a general two or three-dimensional model 
([!]), a corresponding set of coupled Hill's equations may still be obtained by some kind of 
approximation, which can be scaling laws [|I7J] or a variational ansatz ]7[|. 

A new approach involves only imposing the wave-function to have a quadratic dependence 
in the complex phase ip = |^|exp (i J2jk a jkXjXk), j,k = 1,...,3, which is a less restrictive 
ansatz than those used previously. By introducing this trial function into (U), building 
Lagrange's equations for the parameters, and using some transformations, one obtains 

X t + Xt(t)Xi = ^ + Q . (12) 

Here X^,i = 1,2,3 are the three root mean square radius of the solution for each spatial 
direction, i. e. the extensions of the previous X(t) to the three dimensional problem. In 
this case the Aj(t) need not to be equal and in fact experimental problems in BEC involve 
different A values because of asymmetries of the traps. 

Eqs. (|12D are not integrable and form a six-dimensional nonautonomous dynamical 



system for which few things can be said analytically. Nevertheless, the numerical study of 
the approximate model fll2|) exhibits an extended family of resonances which is more or less 
the Cartesian product of those of (PH|), with minor displacements due to the coupling. There 
exist also parameter regions of chaotic and periodic solutions. Numerical simulations of Eq. 
(|I|) confirm the predictions of the simple model (|T2|). 

Nonlinear spectrum and resonances.- Although resonant behavior has been proven in the 
radially symmetric two dimensional case, one could try to make a qualitative explanation 
for the reason that this behavior is also present in non-symmetric problems. Let us look 
for stationary solutions u(r, t) = <f)(r)e~ lflt for the stationary trap, V = V(r), and write the 
nonlinear eigenvalue problem 

/i0= _I A 0+|0| 2 + nr)0. (13) 

Let us assume that we may use these solutions to expand, in a possibly non-orthogonal and 
time dependent way, any solution of the non-stationary problem ([!]). By doing so one finds 
that the energy absorption process in the time dependent potential is ruled by the separation 
between the eigenvalues, /ij — of any two different modes of fllBp. If these differences 
are distributed in a random way, the perturbation is not efficient and does not lead to 
resonances, but to chaos. On the other hand, if the differences may be approximated by 
multiples of a fixed set of frequencies, then an appropriate parametric excitation will induce 
a sustained process of energy gain (and width growth) such as the one we have observed. 

We have studied the spectrum {fi} for the case of an axially symmetric potential V^(r) = 
X r (x 2 + y 2 ) + X z z 2 in three dimensions using a variant of a pseudospectral scheme using the 



harmonic oscillator basis as described in ||18|| . Our results show that the spectrum exhibits 
an ordered structure, with different directions of uniformity that may be excited by the 
parametric perturbation. We have checked the computations in one, two and three spatial 
dimensions without assumptions of symmetry on the solution with similar results despite 
the fact that the spectrum becomes more complex as dimensionality increases. Details of 
these calculations will be published elsewhere. 
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Losses effects. - When losses are included in Eq. (|3]), e.g. by the addition of a new term of 
the type io~u, it is not possible to obtain a set of closed equations for the momenta and thus 
one cannot solve exactly the problem. However, as discussed in the preceding comments 
one still may further restrict the analysis to the parabolic phase approximation to obtain 
the modified equation X + X(t)X = Q(t)/X 3 , where Q(t) is a decreasing function satisfying 
Q — > 1 as t — > oo. At least up to the range of validity of the approximation one obtains 
parametric resonances since the behavior of the singular term does not affect essentially 
the resonance relation. Of course, at the same time one has a decrease in the norm of the 
solution since now dl\/dt = —2al\. 

In conclusion we have studied a parametrically perturbed nonlinear extended system. By 
using the moment technique for the cilindrically symmetric problem we obtain a singular 
Hill equation which is reducible to a linear Hill equation and thus existence of resonances can 
be shown analytically. For the physically relevant case of a periodic perturbation it is shown 
that there exist strong extended resonances for the relevant parameters of the solution even 
when the solution is constrained by conservation laws. It is conjectured using the parabolic 
ansatz, analysis of nonlinear spectrum and numerical simulations that this behavior is also 
present in non-radially symmetric problems or when losses are present. 

V.M.P-G. has been partially supported by CICYT grants PB96-0534 and PB95-0839. P. 
T. is supported by CICYT grant PB95-1203. 
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FIGURES 

FIG. 1. Stability diagram in the e — u) plane. The instability regions corresponding to the first 
higher frequency bands are shaded 
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